/*
 * Equilibrator.h
 *
 *  Created on: 2 Aug 2011
 *      Author: Allan
 */

#ifndef EQUILIBRATOR_H_
#define EQUILIBRATOR_H_

// C++ includes
#include <vector>
#include <string>
using namespace std;

// Eigen includes
#include <Eigen/Core>
using namespace Eigen;

// GeoReact includes
#include "../EquilibriumSolver/EquilibriumSolver.h"

//-----------------------------------------+
//  Model for Condition Specification
//
//  "ChargeBalanced"
//  "Composition::HCO3[-] = 1.0 mmoles"
//  "TotalComposition::Cl[-] = 1.0 moles"
//  "PartialPressure::CO2(g) = 83.2 bar"
//  "Activity::H[+] = 1.0E-7"
//  "pH = 7.2"
//-----------------------------------------+

class Equilibrator 
{
public:
	Equilibrator();
	
	Equilibrator(const vector<string>& allSpecies, const vector<string>& equilibriumSpecies, const vector<string>& equilibriumConditions);
	
	void
	Initialiaze(const vector<string>& allSpecies, const vector<string>& equilibriumSpecies, const vector<string>& equilibriumConditions);
	
	void
	SetTemperaturePressure(double T, double P);
	
	void
	GetState(double&T, double& P, VectorXd& n);
	
	void
	Solve();
	
	void
	Solve(double T, double P, VectorXd& n);
	
	
	/// New Interface
	Equilibrator(const vector<string>& species);
	
	Equilibrator(const Multiphase& multiphase, const ReactionSystem& equilibriumReactions);
	
	void
	Initialiaze(const vector<string>& species);
	
	void
	Initialiaze(const vector<string>& species, const vector<string>& excluded);
	
	void
	AddCondition(const string& type);
	
	void
	AddCondition(const string& type, double value);
	
	void
	AddCondition(const string& type, const string& species, double value);
	
	void
	ClearConditions();
	
	const double
	Get(const string& quantity) const;
	
	const double
	Get(const string& quantity, const string& species) const;
	
	void
	Solve(double T, double P);
	
	void
	PrintState() const;
	
private:
	Multiphase multiphase;
	
	ReactionSystem equilibriumReactions;
	
	vector<EquilibriumCondition> eConditions;
	
	vector<string> equilibriumConditions;
	
	double T, P;
	
	VectorXd n, c, a, g;
};

#endif /* EQUILIBRATOR_H_ */
